Classical theory of universal quantum work distribution in chaotic and disordered non-interacting Fermi systems

We present a universal theory of quantum work statistics in generic disordered non-interacting Fermi systems, displaying a chaotic single-particle spectrum captured by random matrix theory. We consider quantum quenches both within a driven random matrix formalism and in an experimentally accessible microscopic model, describing a two-dimensional disordered quantum dot. By extending Anderson’s orthogonality determinant formula to compute quantum work distribution, we demonstrate that work statistics is non-Gaussian and is characterized by a few dimensionless parameters. At longer times, quantum interference effects become irrelevant and the quantum work distribution is well-described in terms of a purely classical ladder model with a symmetric exclusion process in energy space, while bosonization and mean field methods provide accurate analytical expressions for the work statistics. Our results demonstrate the universality of work distribution in generic chaotic Fermi systems, captured by the analytical predictions of a mean field theory, and can be verified by calorimetric measurements on nanoscale circuits.

where the â i 's stand for fermionic annihilation operators, and the single particle Hamiltonian H(t) incorporates disorder effects and also accounts for the impact of time dependent electrodes. The total fermion number is conserved by Eq. (2), iâ † iâ i = M. We study the universal features of the work statistics for zero temperature quench protocols, starting from the ground state of the initial Hamiltonian, in two steps. First, relying on the observation that the single particle (1) P t (W) ≡ δ W − (Ĥ(t) − E GS (t)) . Deformation-induced motion of energy levels during the quench giving rise to particle-hole excitations. Inset: Quantum quench in a generic disordered dot modelled as a trajectory in the manifold of random matrices. (c) 'Ladder' model for energy absorption: classical diffusion of hard core particles on uniformly spaced energy levels. (d) Benchmarking the work statistics of a realistic quantum dot system against the simplified ladder model for quenches starting from the ground state. Work distribution evaluated from the full microscopic description of the quantum dot, and the results within the classical ladder model collapse onto a universal curve. The values of the numerical parameters are the same as in Fig. 3 below. spectrum of most chaotic systems is captured by random matrix theory (RMT) [44][45][46] , we investigate quenches within random matrix ensembles, such that H(t) follows a trajectory within the manifold of Gaussian random matrices (see the inset of Fig. 1b). Secondly, after establishing the universal properties of work distribution within this random matrix description, we turn to a more realistic microscopic model of a disordered Fermi liquid. We ask the natural question whether the universal features predicted by random matrix theory remain valid for a more realistic, experimentally accessible disordered tight binding Hamiltonian. We find that the answer to this question is affirmative; our results demonstrate the striking universality of the work statistics both within a random matrix framework and for a full microscopic description (see Fig. 1d).
To set up the framework of random matrix theory, we follow the strategy of Refs. 39,47 and 48 , and consider deformations within the space of Gaussian random matrix ensembles, with H 1,2 some independent N × N Gaussian matrices from the orthogonal (GOE), unitary (GUE) or symplectic (GSE) ensembles, inducing the motion and collisions of single-particle energy levels, see Fig. 1b. Here ˙ = v sets the speed of deformations and generates a motion along an 'arc' or ' circle' within the random matrix ensemble, as depicted in the inset of Fig. 1b. Our first goal is to understand universal aspects of the structure and time evolution of the distribution P t (W) within this framework, with the initial state chosen as the ground state of Ĥ (0) . We follow the quantum evolution of the disordered many-body systems, and use a determinant formula presented in Ref. 48 to compute P t (W) . Our first main finding is that the statistics of P t (W) is almost independent of microscopic details as well as the symmetry of the Hamiltonian, once the absorbed energy exceeds sufficiently the one-body energy separation δε ≡ 1/N(ε F ) , characterizing the total density of levels at the Fermi energy ε F , and the time is long enough, t > /δε . To capture work in this long time limit, we construct a simplified classical 'ladder' model which only incorporates Pauli exclusion and level repulsion to the lowest order, in the form of a completely rigid spectrum. However, the 'ladder' model ignores all additional features of the level statistics, including its dependence on the random matrix ensemble parameter β , as well as the interference effects between consecutive level collisions and Landau-Zener transitions, see Fig. 1c. Our simplified 'ladder' model gives a surprisingly accurate description of P t (W) , underlining the high level of universality in the work statistics, and allowing us to derive accurate analytical approximations for P t (W) by means of bosonization and a particle number conserving mean field method. These analytical formulas constitute the second main result of our work. We then turn to our last main objective, and provide further evidence for universality of the work statistics by going beyond the RMT framework. By comparing the 'ladder' model to a realistic 2D quantum dot system, we confirm the striking universality of the work statistics for the full microscopic description, thereby validating our random matrix theory approach (see Fig. 1d).

Results
Quantum mechanical analysis. In this section we briefly summarize the main steps allowing us to numerically evaluate the full distribution of work, P t (W) , for an arbitrary quadratic Hamiltonian, Eq. (2), following Ref. 48 . For a non-interacting Hamiltonian, all information is contained in the time evolution of the single particle wave functions, ϕ m (t) , that evolve according to the Schrödinger equation ( = 1) with initial conditions ϕ m i (0) = δ m i . We solve Eq. (4) using the adiabatic approach, by expanding ϕ m (t) in terms of the instantaneous eigenfunctions η k t of H satisfying H(t) η m t = ε m (t)η m t as We then solve the single particle Schrödinger equation for α m k (t) with initial conditions α m where A kl = −i η k t · ∂ t η l t is the Berry connection. For zero temperature quenches, the characteristic function G t (u) of the work distribution P t (W) can then be expressed by a simple determinant formula 48,49 Here |�(t)� is the full time-evolved many-body wave function, with |�(0� being the ground state of the initial Hamiltonian Ĥ (0), E GS (t) is the ground state energy of the instantaneous Hamiltonian Ĥ (t) , and . . . disord denotes the averaging over disorder. The matrix g t (u) contains information on overlaps and the instantaneous single particle energies ε k (t) at time t, We compute g t (u) numerically, average over disorder or the random matrix ensemble, . . . disord , and determine the final distribution by performing a Fourier transformation.
The average level spacing δε and its inverse provide natural energy and time scales, and allow us to introduce the dimensionless work and time, w ≡ W/δε and t ≡ t δε , respectively. As shown in Fig. 1b, deformations of www.nature.com/scientificreports/ the Hamiltonian lead to a continuous motion of single particle levels, and thereby induce collisions and transitions between them. At longer times, these collisions and Landau-Zener transitions 38,40,47,50 give rise to a diffusive broadening of the Fermi surface, more precisely, of the Fermi edge in energy space, since momentum is not well-defined in the presence of disorder. Namely, the mean occupation of the kth energy level is given by f k (t) = (1 − t k (t))/2 with t k (t) = erf �k/ 4 Dt where �k = k − M − 1/2 is measured from the Fermi-level, and D is the dimensionless energy diffusion constant. After a short time perturbative ∼ t 2 scaling, the average work is found to increase as �w� = Dt (see Refs. 37,38,40,48 and Methods). Let us note that quantum effects such as dynamical localization yielding time dependent corrections to the diffusion constant studied for periodic driving in the important early works 41-43 become only relevant on a time scale much larger than the period 39 and are thus negligible on our time scales. Below we turn to the full distribution of the dimensionless work w in different settings, both numerically and through various analytical approximations.
Quenches in random matrix ensembles. In this section we analyze the works statistics within the framework of random matrix theory for quenches of the form of Eq. (3). The disorder average appearing in Eq. (6) becomes an average over two independent random matrices, H 1,2 . The distribution P˜t (w) can be disentangled into an adiabatic and a regular part, Random matrix theory implies that -apart from the symmetry of the Hamiltonian -the statistics of the evolution of the eigenvalues, sketched in Fig. 1b, is completely characterized by the velocity with which levels deform, i.e., the frequency of avoided level crossings. Indeed, the average distance of level crossings 47,50,51 , , and the time scale 1/δε define a natural 'velocity' in parameter space, v c ≡ �� �δε , which we can use to introduce the dimensionless velocity, ṽ ≡˙ /(�� �δε) . For N × N random matrices, δε ∼ 1/N , and �� � ∼ 1/ √ N , therefore v c ∼ 1/N 3/2 . The dimensionless velocity characterizes microscopic processes. For ṽ ≪ 1 the motion is almost adiabatic, and small probability Landau-Zener transitions dominate. For ṽ ≫ 1 , on the other hand, transitions between remote levels generate energy absorption.
From our random matrix considerations it follows that the distribution P˜t (w) can only depend on t , ṽ , and, in case of finite temperature initial states, on the dimensionless initial temperature, T ≡ T/δε . Similarly, the diffusion constant D is a universal function of ṽ , which scales as D ∼ṽ 2 for large velocities, while for ṽ < 1 nearest neighbor transitions dominate and yield D ∼ṽ (β/2+1) , with β = 1, 2 and 4 characterizing the orthogonal, unitary, and symplectic ensembles, respectively (see the Methods for more details).
We show our random matrix simulation results in Fig. 2. We checked that the work distribution is not sensitive to changes in the number of energy levels N and the number of particles M, as long as the particle-hole excitations contributing to the work are created far from the edges of the band (see also Ref. 48 .) For small work, w 10 , www.nature.com/scientificreports/ the statistics depends on β as well as on ṽ, and P reg (w;t ) displays peaks and minima associated with level repulsion, clearly reflecting the symmetry of the underlying Hamiltonian (see Fig. 2a). These features become more pronounced for larger β due to the stronger level repulsion. For larger works, �w� max {ṽ 2 , 1} , however, one enters a diffusion dominated regime, where symmetry related and microscopic features become less important, and a universal distribution displayed in Fig. 2b emerges. The observed distribution is clearly non-Gaussian, and characterizes work statistics in generic non-interacting fermion systems.
Ladder model. The agreement between the three universality classes is suggestive that quantum interference effects do not play an important role in this diffusion-dominated regime. We can therefore attempt and construct a classical 'ladder' model, consisting of uniformly placed classical energy levels at a distance δε from each other (see Fig. 1c), occupied by hard-core particles in line with Fermi statistics. The energy of a many-body state is then given by E = k n k ε k with n k ∈ {0, 1} the occupation numbers, and k n k = M the total number of particles. The evenly placed levels (9) mimic level repulsion and level rigidity in chaotic systems. As a final component, perturbationinduced random Landau-Zener transitions are modeled by nearest neighbor hopping transitions and a symmetrical exclusion process (SEP) in energy space. The SEP is a classical Markovian stochastic system of hard-core particles trying to jump at a constant rate to the left or to the right with probability 1/2 along a 1D lattice with the constraint that no site can be occupied by more than one particle (exclusion property). This simple model captures the diffusive broadening of the Fermi surface (see Methods) and, in addition to level repulsion, it also incorporates Fermi statistics and particle number conservation. In our SEP simulations, the jump rate was set by the diffusion constant D such that for a given quench time the same final occupation profile and average work is obtained as in the quantum simulation. As can be seen in Figs. 2 and 3, this classical stochastic model gives a surprisingly accurate description of the work statistics for large enough average work, independently of the velocity. Moreover, with certain assumptions, the 'ladder' model can be used to compute P ad (t ) and P reg (w;t ) analytically for a T = 0 temperature initial state, without performing the actual Monte Carlo simulations, using either bosonization or a more accurate mean field approach. It is, however, crucial to treat particle number conservation with care.
Bosonization. Bosonization offers a simple method to treat particle number conservation in the 'ladder' model. Introducing fermion operators for each level, we can express the total energy as H = k (ε k − ε F ) : c † k c k : with ε F = δε (M + 1/2) the Fermi energy and : ... : referring to normal ordering with respect to the Fermi sea. Following Ref. 52 q ′ ] = δ q,q ′ , and rewrite the Hamiltonian in terms of these as with N = k c † k c k − M the normal ordered fermion number. Clearly, the fermion number does not change for the closed system studied here so the second term in Eq. (10) does not give a contribution. We can obtain an approximate expression for P˜t (w) by assuming that the final state is thermal with an effective boson temperature T eff = √ 6�w�/π , chosen to yield the appropriate average energy, � q>0 q b † q b q � ≡ �w� . In the large w limit, we then obtain (for details see Methods),  www.nature.com/scientificreports/ where I 1 is the modified Bessel function of the first kind. Since T eff ∼ Dt , the prefactor decays as ∼ e −C √ Dt , corresponding to a stretched exponential decay of adiabatic processes, as confirmed by our quantum mechanical simulations 48 .
Mean field theory. The bosonization approach yields a good account of the overall structure of P˜t (w) , but with certain limitations (see Fig. 2b). In particular, the assumption of a thermal final state is not quite correct. The occupation of the single particle levels after the time evolution is not described by the Fermi function but has a diffusive structure, as stated earlier. A more accurate expression can be obtained for P˜t (w) in a simple, particle number conserving mean field approach, where instead of assuming thermalization, we rely on the diffusive nature of energy absorption, and assume that each fermion level k is occupied with probability f k = (1 − erf �k/ 4 Dt )/2 with �k = k − M − 1/2 measured from the Fermi-level, corresponding to a diffusive broadening of the Fermi surface. To enforce the constraint, k n k = M , we use an integral representation over an auxiliary variable. A saddle point procedure in this latter then yields accurate expressions for P ad (t ) as well as for P reg (w;t ).
The mean field probability distribution, P MF t (w) , is similar in structure to Eq. (11), but contains additional correction terms (see Methods for details), with c w ≈ 1.35 and As shown in Fig. 2b, the mean field expressions above yield an accurate description of work in the diffusive regime. Similar to the bosonization result, Eq. (11), P MF t (w) is non-Gaussian and, by construction, depends parametrically only on w . The probability of adiabatic processes also falls off as a stretched exponential, but the prefactor c w is more accurate than the one obtained by the simple bosonization theory ( π/ √ 6 ≈ 1.28 ) (see Ref. 48 and Methods).

Validation by microscopic models and experimental setup.
The universal features of the work statistics discussed so far have been established within the framework of random matrix theory. It is a natural question to that extent do these result remain valid outside this framework, in realistic disordered systems. To answer this question and to confirm our predictions as well as to validate the results of our random matrix approach, we propose to study a 2-dimensional quantum dot, and squeeze the electron gas confined there by applying time dependent external gate voltages (see Fig. 1a). This system can be realized experimentally 53,54 .
We model the quantum dot by a disordered tight binding Hamiltonian on a square lattice, where the sum over δ runs over the two orthogonal lattice vectors. The first term accounts for the kinetic energy of the electrons, while ε r are random onsite energies drawn from a Gaussian distribution of variance σ , and are responsible for electron scattering and disorder. We emphasize that we focus on the delocalized regime of the model, where the Anderson localization length is much larger than the system size. The potential describes the parabolic confinement, generated by external gate electrodes. The second term in V (r, t) describes a compression (decompression) of the electron gas in the x direction with a simultaneous decompression (compression) along the y direction. We vary linearly in time between − f and f to induce deformations and generate dissipation. A numerical investigation of the single particle spectrum of Eq. (14) reveals that, although some deviations are clearly present, the spectrum of Eq. (14) is reasonably described in terms of GOE for each value of (see Methods).
Similarly to the RMT simulations, we generate work by varying lambda at a constant pace, ˙ = v, corresponding to the velocity of potential deformations, and compute work statistics P˜t (w) according to the formula used for the random matrix models, Eq. (6). In order to get universal work distribution, we introduce the same velocity and time units, however, their dependence on the number of lattice sites, on-site disorder strength, and potential strength is quite different from that in the RMT model. In particular, both δε and the typical distance between avoided crossings strongly depend both on the potential strength, system size, on-site disorder strength, and the position in the energy spectrum which we, therefore, computed numerically, averaging over ∼ 5 × 10 3 disorder realizations. Then setting the average work in these units to �w� = 10, we computed the regular part of the work statistics, P reg (w;t) , for velocity ṽ = 0.4 , presented in Fig. 3. www.nature.com/scientificreports/ The results show striking agreement with random matrix theory as well as with the 'ladder' model, and thereby provide further evidence for universality.
An alternative experimental platform to study quantum work statistics is offered by ultracold atoms 17 . For a forward-backward protocol P ad is essentially the ground state fidelity, which has been measured in Ref. 55 by preparing two identical copies of a quantum system, and measuring their overlap. This method could be used to verify the predicted stretched exponential behavior of P ad in disordered fermion systems.

Discussion
We demonstrated universal behavior in the quantum thermodynamics of a chaotic Fermi liquid, displaying a level statistics captured by random matrix theory, by studying the work statistics for quantum quenches in disordered non-interacting fermionic systems. We note that according to Fermi liquid theory, interaction effects are expected to remain irrelevant at single-particle energies ε such that (ε − ε F ) 2 /ε F is smaller than the typical level spacing.
First we focused on quenches formulated within the framework of random matrix theory. Strikingly, we found that for large enough average work, the distribution is independent of the random matrix ensemble and is very well captured by a universal classical stochastic model describing diffusion in energy space. This simplified 'ladder' model only relies on Pauli exclusion principle and on the rigidity of the random matrix spectrum, but neglects all further microscopic details, including the symmetry class of the Hamiltonian. This high level of simplification allowed us to derive approximate analytical expressions both via bosonization and mean field theory. Interestingly, the bosonization result in Eq. (11) also emerged in the context of work statistics in Luttinger liquids after an interaction quench 26 . While the bosonization approach performs more poorly in comparison with the mean field treatment (cf. Fig. 2b), it relies on the additional approximation that the final state is thermal. In contrast, the mean field approach incorporates the more accurate diffusive occupation profile of the final state, and yields a striking agreement with the results of full quantum simulations, thereby providing an excellent analytical approximation for the universal work statistics. Finally, we turned to a more realistic, experimentally accessible microscopic model, and studied the work distribution of a 2D disordered quantum dot subject to timedependent gate voltages. We demonstrated the striking agreement between the work statistics of this microscopic model and our random matrix theory predictions. These results show that our universal description holds beyond the framework of random matrix theory, and accurately captures the energy absorption in realistic models. Our results could be tested experimentally by studying squeezed disordered quantum dots.

Bosonization approach
In this approach, we consider an equilibrium fermionic system with uniformly spaced one-particle energy levels. In the framework of bosonization, the fermionic particle-hole excitations with respect to the ground state are represented as bosonic states. We assign thermal Boltzmann weights e −β eff qδε to these states, where β eff is an effective inverse temperature while q δε with q = 1, 2, . . . measures the energy of the particle-hole excitation.
Since these excitations are bosonic, for each q we can have n q = 0, 1, 2, . . . arbitrarily many bosonic excitations with energy qn q δε . In the characteristic function each of them carries a contribution of e iũqn q , so we have where N = q>0 1 − e −q(β eff δε) −1 so that G eff (0, T eff ) = 1 . Exponentiating Eq. (16) and taking the continuum limit q>0 → ∞ 0 dx we get: The Fourier transform can be performed exactly, which together with the normalisation factor leads to the analytical result expressed in terms of the dimensionless effective temperature Mean field approach. In this section we provide some details about the mean field theory calculations and the resulting analytic expressions.
6 w e −w/ T eff I 1 π 2 3 w + δ(w) . www.nature.com/scientificreports/ Probability of adiabaticity. Within the mean field approach, the probability of each many-body configuration takes the form of the product of independent Bernoulli weights of M occupied and N − M empty sites. In order to simplify calculations and without any loss of generality we consider the case of M = N/2: where the particle number conservation is taken into account by the Kronecker-delta for which we used a standard integral representation. The Bernoulli weights are where f k (t) = (1 − t k (t))/2 with t k (t) = erf �k/ 4 Dt and �k = k − M − 1/2 is measured from the Fermilevel. Finally, the time-dependent normalization factor is the sum of all possible many-body probabilities: Writing the above expression as the exponential of its logarithm, approximating the resulting sum by an integral and performing a saddle point approximation around = 0, we obtain for large enough values of �w� = Dt ≫ 1: The probability of adiabaticity then reads with C ≈ 1.35.

Variance of work.
For �w� ≫ 1 , we approximate the variance of the work by neglecting the fluctuations of the energy levels 56 but incorporating the fluctuations of the occupation numbers. We thus measure the energies from the Fermi-level as ε k (t) → �k δε with �k = k − M − 1/2. For a given realization of H(t) , this leads to the estimate where . . . denotes quantum average. Separating the diagonal terms, the RM average δw 2 (t) RM can be written as where δn k,t ≡n k,t − �n k,t � is the deviation of the occupation number from the mean value. As the n k,t behave as binary random variables, the averages in the first term are given by ��δn 2 The correlators in this equation can be expressed in terms of the amplitudes α m k (t) as �δn k,t δn k ′ ,t � = − The negativity of this correction implies that the level occupations are anticorrelated, as follows from particle number conservation. Neglecting this correction for the moment and replacing sums by integrals, we arrive at the estimate yielding �δw 2 (t)� ∼ �w� 3/2 . We thus recovered the observed behavior, however, the prefactor turns out to be incorrect. The correct prefactor 57 can be obtained by a more careful mean field calculation that takes into account the occupation number correlations as well which obey that same scaling ∼t 3/2 .
Distribution of work. The characteristic function of the distribution of work can be expressed as www.nature.com/scientificreports/ where we introduced the scaled variable ũ = u δε and the notation h k (u, t) = Here . . . MF denotes averaging over the mean field many-body probabilities and N t a modified normalization constant. As numerics revealed, for large enough injected works �w� ≫ 1 neglecting particle number conservation does not introduce big errors provided we subtract the pure particle and pure hole excitations with respect to the ground state: where the first term is the = 0 saddle-point solution of the integral expression, while the second part substracts the contributions coming from the pure particle-hole excitations. Here the integrals can be approximated as yielding with c w = 3 √ 2π 5 chosen such that the characteristic function correctly reproduces the first two cumulants of work in the saddle point solution. Now, analogously to the Bosonization approach, both terms can be Fourier transformed exactly leading to the approximate analytic expression Energy space diffusion. In this section we demonstrate that the energy level occupations exhibit a diffusive profile, meaning that particle-hole excitations happen dominantly in a window growing as ∼ �w� 1/2 , for all the random matrix ensembles as well as for the "ladder model" and the disordered quantum dot. The left panel of Fig. 4 shows that for large enough average work the mean level occupation for of all three RMT ensembles (GOE, GUE, GSE) follows a single universal curve identical to those of the quantum dot model up to high precision and it is also perfectly described by the ladder model. Numerical calculations were made for ∼ 5 × 10 3 disorder realizations both for random matrix theory and the disordered quantum dot, for N = 40, 28, 40 for the three ensembles, respectively and for parameters L = 38, σ = 1.75J, α = 75J and with 427 particles in the case of the quantum dot.
The right panel of Fig. 4 shows the velocity dependence of the diffusion constant, D β (ṽ) for the three ensembles and the quantum dot model. We averaged over ∼ 5 × 10 3 simulations, yielding smooth enough timeevolutions of average work to extract the diffusion constants. Parameters were chosen such that we avoid finite size effects and be in the diffusion regime. The rate of energy absorbed by the system exhibits an anomalous frequency dependence for slow quenches, D β (ṽ 1) ∼ṽ β/2+1 , while for fast processes becomes independent of the underlying symmetry class and grows quadratically, as it should in the case of a metal, D β (ṽ ≫ 1) ∼ṽ 2 . The diffusion constant for the quantum dot shows the same power-law behavior as the GOE ensemble, albeit with a slightly smaller prefactor. Finally, we compare the level spacing distribution of the GOE ensemble and the disordered quantum dot. As shown in Fig. 5, the distribution of the distance of neighboring levels are well described by the analytical RMT result given by the Wigner surmise. Similar observations hold for the statistics of the the Landau-Zener parameters at the avoided level crossings in comparison with the RMT results of Ref. 37 .

Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.  Figure 4. Energy space diffusion. Left: Average occupations of instantaneous single particle eigenstates for the three RMT ensembles (GOE, GUE, GSE) and the quantum dot model compared to the classically obtained results within the ladder model. All the five curves collapse onto a single universal, diffusively broadening profile given by 1 − erf �k/ √ 4�w� /2 . Right: Velocity dependence of the diffusion constant. For slow quenches it has an anomalous power-law behavior, D(ṽ 1) ∼ṽ β/2+1 , while for fast quenches it grows quadratically and with the same prefactor for the RMT ensembles. The quantum dot displays a similar behavior as the GOE ensemble in the two limiting cases, with a slightly different prefactor. In the left panel, we used N = 40 energy levels for the GOE and GSE, N = 28 for the GUE ensemble, and N = 56 for the ladder model simulations, and half filling in all cases. For the quantum dot model, the same parameters were used as in Fig. 3. We used various values of RMT and quantum dot parameters to obtain the results shown in the right panel. , obtained by Wigner's surmise describing the GOE case. For the numerical calculations we averaged over ∼ 5 × 10 4 disorder realizations which proved to yield smooth enough curves.